#include "element.h"
#include<QDebug>
#include <QVector>

using namespace arma;


element::element()
{
    lambda = 1.5;
    mi = 1.0;

    setGaussPoint();
}

void element:: setNumeration(int num)
{
    numeration = num;
}

int element:: getNumeration()
{
    return numeration;
}

void element:: setNodes(node* n1, node* n2, node* n3 , node* n4)
{
    nos.clear();
    nos.append(n1);
    nos.append(n2);
    nos.append(n3);
    nos.append(n4);
}

node* element::getNode(int i)
{
    return nos[i-1];
}


void element::setLocalMatrix()
{
    localMatrix=zeros(8,8);
    for (int a=0 ; a<4 ; a++)
    {
        for (int b=0 ; b<4 ; b++)
        {
            for (int i=0 ; i<2 ; i++)
            {
                for (int j=0 ; j<2 ; j++)
                {
                    for (int r=0 ; r<2 ; r++)
                    {
                        for (int s=0 ; s<2 ;s++)
                        {
                            localMatrix(2*(a)+r,2*(b)+s) += getTensorCComponent(i,r,j,s)*getTensorQabijComponent(a,b,i,j);
                        }
                    }
                }
            }
        }
    }
//    localMatrix = ones(8,8);
}

void element::setLocalForce()
{
    localForce = zeros(8,1);
    mat aux;
    for(int a=0 ; a<4 ; a++) //f contribution
    {
        for(int r=0 ; r<2 ; r++)
        {
            for(int b=0 ; b<4 ; b++)
            {
                localForce(2*a+r,0) += f(r,b)*getTensorQabComponent(a,b);
            }
        }
    }


    for(int r = 0 ; r<2 ; r++)//prescribed contribution
    {
        for(int b = 0 ; b<4 ; b++)
        {
            aux = nos[b]->getBoundaryCondition();
            if(aux(0,r)==1)
            {
                for (int m =0 ; m<8 ; m++)
                {
                    localForce(m,0) -= localMatrix(m,2*b+r)*aux(1,r);
                }
            }
        }
    }
}

mat element::getLocalMatrix()
{
    return localMatrix;
}

vec element::getLocalForce()
{
    return localForce;
}

float element::f(int r, int b) // has to be modified for each problem
{
    float x = nos[b]->getPositionX();
    float y = nos[b]->getPositionY();
    if(r==0)
    {
        return 2*(cos(datum::pi*x)*cos(datum::pi*y));
    }
    else if(r==1)
    {
        return 2*(sin(datum::pi*x)*sin(datum::pi*y));
    }
    else
    {
        return 0.0;
    }
}

float element::getTensorCComponent(int i, int r, int j ,int s)
{
    if((i==j) && (j==r) && (r==s))
    {
        return lambda + 2*mi;
    }
    else if ((i==r) && (j==s))
    {
        return lambda;
    }
    else if(((i==j) && (r==s)) || ((i==s)&&(r==j)))
    {
        return mi;
    }
    else
    {
        return 0;
    }
}


void element::setJacobiano(float ksi, float eta)
{
    jacobiano = zeros(2,2);
    setPhi(ksi,eta);
    setDphi(ksi,eta);

    jacobiano(0,0) = getNode(1)->getPositionX()*dPhi(0,0) + getNode(2)->getPositionX()*dPhi(1,0) + getNode(3)->getPositionX()*dPhi(2,0) + getNode(4)->getPositionX()*dPhi(3,0);
    jacobiano(0,1) = getNode(1)->getPositionX()*dPhi(0,1) + getNode(2)->getPositionX()*dPhi(1,1) + getNode(3)->getPositionX()*dPhi(2,1) + getNode(4)->getPositionX()*dPhi(3,1);
    jacobiano(1,0) = getNode(1)->getPositionY()*dPhi(0,0) + getNode(2)->getPositionY()*dPhi(1,0) + getNode(3)->getPositionY()*dPhi(2,0) + getNode(4)->getPositionY()*dPhi(3,0);
    jacobiano(1,1) = getNode(1)->getPositionY()*dPhi(0,1) + getNode(2)->getPositionY()*dPhi(1,1) + getNode(3)->getPositionY()*dPhi(2,1) + getNode(4)->getPositionY()*dPhi(3,1);
}



mat element::getJacobiano()
{
    return jacobiano;
}

void element::setPhi(float ksi, float eta)
{
    phi = zeros(4,1);
    phi(0,0) = (1-ksi)*(1-eta)/4;
    phi(1,0) = (1+ksi)*(1-eta)/4;
    phi(2,0) = (1+ksi)*(1+eta)/4;
    phi(3,0) = (1-ksi)*(1+eta)/4;
}


void element::setDphi(float ksi, float eta)
{
    dPhi = zeros(4,2);
    dPhi(0,0) = -(1-eta)/4;
    dPhi(1,0) =  (1-eta)/4;
    dPhi(2,0) =  (1+eta)/4;
    dPhi(3,0) = -(1+eta)/4;
    dPhi(0,1) = -(1-ksi)/4;
    dPhi(1,1) = -(1+ksi)/4;
    dPhi(2,1) =  (1+ksi)/4;
    dPhi(3,1) =  (1-ksi)/4;
}

mat element::getPhi()
{
    return phi;
}


mat element::getDphi()
{
    return dPhi;
}


void element::setGaussPoint()
{
    gaussPoint = zeros(4,2);

    float a = sqrt(3)/3;

    gaussPoint(0,0) = -a;
    gaussPoint(0,1) = -a;
    gaussPoint(1,0) =  a;
    gaussPoint(1,1) = -a;
    gaussPoint(2,0) =  a;
    gaussPoint(2,1) =  a;
    gaussPoint(3,0) = -a;
    gaussPoint(3,1) =  a;
}

float element:: getTensorQabijComponent(int a,int b ,int i ,int j)
{
    float value = 0.0;
    mat inverseJacobTran;


    for(int t=0; t< 4; t++)
    {
        setJacobiano(gaussPoint(t,0),gaussPoint(t,1));
        inverseJacobTran  = jacobiano.t().i();
        value += (dPhi(a,0)*inverseJacobTran(i,0)+dPhi(a,1)*inverseJacobTran(i,1))*(dPhi(b,0)*inverseJacobTran(j,0)+dPhi(b,1)*inverseJacobTran(j,1))*det(jacobiano);
    }

    return value;
}

float element::getTensorQabComponent(int a, int b)
{
    float value = 0.0;

    for (int t=0; t<4 ; t++)
    {
        setJacobiano(gaussPoint(t,0),gaussPoint(t,1));
        value += phi(a,0)*phi(b,0)*det(jacobiano);
    }
    return value;
}

QDebug operator<<(QDebug x,element e)
{
    x.nospace() << "elemento[" << e.getNumeration() << "]: "<< *e.getNode(1)<< " " << *e.getNode(2) << " " << *e.getNode(3) << " " << *e.getNode(4);
    return x.space();
}



